Setup relevant libraries and directory for saving figs and data

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(ggplot2)
library(tidyr)
library(zoo)

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

devtools::load_all('..')

## ℹ Loading hector

data_dir <- 'data/'
ini_dir <- 'data/pic_hector_ini/climate/'

GCP data for NBP

Read in and understand GCP data

This is true NBP: atmosphere to land, positive NBP = a sink

# load the file passed from Dawn
gcp_data_dawn <- read.csv(paste0(data_dir, "nbp_gcp.csv"))

# Load the file pulled directly from GCP:
read.csv(paste0(data_dir,'Global_Carbon_Budget_2022v1p0_historical_co2_tab.csv'), stringsAsFactors = F) %>%
  select(year=Year, e_luc=land.use.change.emissions, s_land = land.sink, Units) %>%
  tidyr::replace_na(list(e_luc=0)) %>%
  # If we want to sum e_luc + s_land, uncomment from the `gather` to `ungroup` 
  # calls here and comment out the `mutate`:
  # gather(reporting_id, value, -year, -Units) %>%
  # group_by(year, Units) %>%
  # summarize(nbp = -sum(value, na.rm=T)) %>%
  # ungroup %>%
  # Following slide 56 of GCP_data/papers/GCP_CarbonBudget_2023_slides_v1.0-2-Alissa-Haward.pdf
  mutate(nbp = (s_land - e_luc)) ->
  gcp_data

Check if GCP data and Dawn’s agree

Check if nbp = (s_land-e_luc) from the raw GCP download agrees with file from Dawn. Note Dawn’s file did -(s_land-e_luc) to have things in the land to atmosphere direction instead, and compare with offline model outputs directly. so have to do minus a minus

print("Difference to Dawn's file:")

## [1] "Difference to Dawn's file:"

print(max(abs(gcp_data$nbp - -gcp_data_dawn$nbp)))

## [1] 0.01

Proceed with GCP data

Label data and convert units

# Label data and convert units
gcp_data$scenario <- "1. Global Carbon Project"
gcp_data$nbp_raw <- gcp_data$nbp*1000
gcp_data$nbp <- zoo::rollmean(gcp_data$nbp_raw,k=10, align='right', fill= NA)

gcp_data %>% 
    filter(year >=1970)

##    year e_luc s_land  Units  nbp                 scenario nbp_raw
## 1  1970  1.26   0.71 GtC/yr -228 1. Global Carbon Project    -550
## 2  1971  1.25   2.45 GtC/yr  -25 1. Global Carbon Project    1200
## 3  1972  1.25   1.25 GtC/yr    4 1. Global Carbon Project       0
## 4  1973  1.21   1.89 GtC/yr  139 1. Global Carbon Project     680
## 5  1974  1.18   4.19 GtC/yr  420 1. Global Carbon Project    3010
## 6  1975  1.19   2.65 GtC/yr  670 1. Global Carbon Project    1460
## 7  1976  1.17   2.99 GtC/yr  844 1. Global Carbon Project    1820
## 8  1977  1.20   1.51 GtC/yr  841 1. Global Carbon Project     310
## 9  1978  1.15   2.78 GtC/yr  900 1. Global Carbon Project    1630
## 10 1979  1.11   1.42 GtC/yr  987 1. Global Carbon Project     310
## 11 1980  1.11   0.41 GtC/yr  972 1. Global Carbon Project    -700
## 12 1981  1.20   2.22 GtC/yr  954 1. Global Carbon Project    1020
## 13 1982  1.20   1.67 GtC/yr 1001 1. Global Carbon Project     470
## 14 1983  1.31   0.28 GtC/yr  830 1. Global Carbon Project   -1030
## 15 1984  1.45   2.89 GtC/yr  673 1. Global Carbon Project    1440
## 16 1985  1.39   2.77 GtC/yr  665 1. Global Carbon Project    1380
## 17 1986  1.40   2.31 GtC/yr  574 1. Global Carbon Project     910
## 18 1987  1.40   0.35 GtC/yr  438 1. Global Carbon Project   -1050
## 19 1988  1.38   2.06 GtC/yr  343 1. Global Carbon Project     680
## 20 1989  1.33   3.57 GtC/yr  536 1. Global Carbon Project    2240
## 21 1990  1.33   2.33 GtC/yr  706 1. Global Carbon Project    1000
## 22 1991  1.30   2.25 GtC/yr  699 1. Global Carbon Project     950
## 23 1992  1.33   2.26 GtC/yr  745 1. Global Carbon Project     930
## 24 1993  1.33   2.86 GtC/yr 1001 1. Global Carbon Project    1530
## 25 1994  1.45   1.47 GtC/yr  859 1. Global Carbon Project      20
## 26 1995  1.44   1.77 GtC/yr  754 1. Global Carbon Project     330
## 27 1996  1.44   3.40 GtC/yr  859 1. Global Carbon Project    1960
## 28 1997  1.94   3.19 GtC/yr 1089 1. Global Carbon Project    1250
## 29 1998  1.50   1.51 GtC/yr 1022 1. Global Carbon Project      10
## 30 1999  1.50   3.59 GtC/yr 1007 1. Global Carbon Project    2090
## 31 2000  1.40   3.75 GtC/yr 1142 1. Global Carbon Project    2350
## 32 2001  1.29   2.38 GtC/yr 1156 1. Global Carbon Project    1090
## 33 2002  1.41   0.99 GtC/yr 1021 1. Global Carbon Project    -420
## 34 2003  1.54   2.31 GtC/yr  945 1. Global Carbon Project     770
## 35 2004  1.46   3.47 GtC/yr 1144 1. Global Carbon Project    2010
## 36 2005  1.29   1.88 GtC/yr 1170 1. Global Carbon Project     590
## 37 2006  1.38   3.11 GtC/yr 1147 1. Global Carbon Project    1730
## 38 2007  1.21   2.71 GtC/yr 1172 1. Global Carbon Project    1500
## 39 2008  1.27   3.63 GtC/yr 1407 1. Global Carbon Project    2360
## 40 2009  1.37   2.99 GtC/yr 1360 1. Global Carbon Project    1620
## 41 2010  1.32   3.32 GtC/yr 1325 1. Global Carbon Project    2000
## 42 2011  1.35   4.11 GtC/yr 1492 1. Global Carbon Project    2760
## 43 2012  1.32   2.32 GtC/yr 1634 1. Global Carbon Project    1000
## 44 2013  1.26   3.56 GtC/yr 1787 1. Global Carbon Project    2300
## 45 2014  1.34   3.66 GtC/yr 1818 1. Global Carbon Project    2320
## 46 2015  1.47   2.01 GtC/yr 1813 1. Global Carbon Project     540
## 47 2016  1.24   2.69 GtC/yr 1785 1. Global Carbon Project    1450
## 48 2017  1.18   3.56 GtC/yr 1873 1. Global Carbon Project    2380
## 49 2018  1.14   3.65 GtC/yr 1888 1. Global Carbon Project    2510
## 50 2019  1.24   3.04 GtC/yr 1906 1. Global Carbon Project    1800
## 51 2020  1.11   3.11 GtC/yr 1906 1. Global Carbon Project    2000
## 52 2021  1.08   3.45 GtC/yr 1867 1. Global Carbon Project    2370

compare hector NBP with GCP

#comparison with GCP

# package scenario
rcp <- "ssp245"
scenario_file <- paste0("input/hector_",rcp,".ini")
ini_file245 <- system.file(scenario_file, package="hector") 
#### QUESTION - is this right hector package if did devtools::load_all
#### ANSWER: yes, because added different ini to inst/input directory and after
####         devtools::load_all(..) that new ini is then accessible.

core245 <- hector::newcore(ini_file245, suppresslogging = FALSE, name='245')
hector::run(core245,runtodate = 2100)

## Hector core: 245
## Start date:  1745
## End date:    2300
## Current date:    2100
## Input file:  /Users/snyd535/Documents/GitHub/hector/inst/input/hector_ssp245.ini

out <- hector::fetchvars(core245, dates=1746:2100, vars=c(hector::NBP()))

LUC emissions received by Hector

gcam_emiss <- read.csv('/Users/snyd535/Documents/GitHub/hector/inst/input/gcam_emissions.csv', skip =4)

ggplot(gcam_emiss) +
    geom_line(mapping = aes(x=Date, y=luc_emissions))

library(rgcam)

get_gcam_luc_emissR <- function(db_name="database_basexdb", gcam_dir, scenario,
                                read_from_file=FALSE,
                                filename="data/gcam_land_alloc.csv"){
  
  if (read_from_file) {
    gcam_land_alloc <- read.csv2(file=filename,header=TRUE)
  }
  else {
    base_conn <- rgcam::localDBConn(gcam_dir, db_name)
    
    land_alloc_query <- '<query title="LUC emissions by region">
                <axis1 name="LandLeaf">LandLeaf</axis1>
                <axis2 name="Year">land-use-change-emission[@year]</axis2>
                <xPath buildList="true" dataName="land-use-change-emission" group="false" sumAll="true">/LandNode[@name=\'root\' or @type=\'LandNode\' (:collapse:)]//
                land-use-change-emission[@year&gt;1970]/text()</xPath>
                <comments/>
            </query>'
    
    new.proj <- rgcam::addSingleQuery(base_conn, "new.proj", "LUC_R", land_alloc_query, c(scenario), clobber=TRUE)
    
    gcam_land_alloc <- rgcam::getQuery(new.proj, "LUC_R")
   # write.csv(gcam_land_alloc, file=filename, row.names = FALSE)
  }
  
  return(gcam_land_alloc)
}


luc_R<- get_gcam_luc_emissR(db_name = 'database_basexdbGCAM',
                                 gcam_dir = 'data/pic_hector_ini',
                               scenario = 'GCAM')

## Database scenarios:  GCAM

ggplot(luc_R %>% filter(landleaf == 'Pasture_MissouriR', region == 'USA', year <=2015)) +
    geom_line(mapping = aes(x=year, y=value, color = scenario)) +
    ylab('MtC/yr')+
    ggtitle('Pasture_MissouriR GCAM output LUC emissions')

ggplot(luc_R %>% 
           group_by(Units, scenario, year) %>%
           summarize(value = sum( value)) %>%
           ungroup %>% 
           filter(year<=2015)) +
    geom_point(mapping = aes(x=year, y=value, color = scenario)) +
    ylab('MtC/yr')+
    ggtitle('Global aggregate GCAM output LUC emissions') +
    geom_line(data = gcam_emiss %>% filter(Date >= 1975, Date <= 2015),
              mapping = aes(x=Date, y=luc_emissions*1e3), color = 'black')

## `summarise()` has grouped output by 'Units', 'scenario'. You can override using
## the `.groups` argument.

Plot NBP data

gcp_data %>%
  select(year, scenario, nbp, Units) %>%
    select(year, value=nbp, units=Units) %>%
    mutate(variable = 'NBP',
           scenario ="GCP",
           value = 0.001*value,
           units = 'PgC/year'
           ) %>% 
  bind_rows(out) ->
  out2


ggplot(out2 %>% filter(year <= 2021)) +
    geom_line(mapping = aes(x=year,y=value, colour=scenario)) +
    geom_line(data= (gcp_data %>%
                         select(year, scenario, nbp_raw, Units) %>%
                         select(year, value=nbp_raw, units=Units) %>%
                         mutate(variable = 'NBP_raw',
                                scenario ="GCP_raw",
                                value = 0.001*value,
                                units = 'PgC/year'
                                )),
              mapping = aes(x=year, y=value), color = 'grey', alpha = 0.5) +
    ggtitle('gray is raw GCP')

## Warning: Removed 9 rows containing missing values or values outside the scale range
## (`geom_line()`).

rm(out)
rm(out2)
rm(gcp_data_dawn)
rm(luc_R)

Test Hector NBP constraints.

The 0-level test for whether the NBP constraint works is if you can do a Hector run, extract the NBP from that run, use that NBP as a constraint for a new Hector run, and get the same results.

This is trickier than it sounds, due to uninformative errors triggered in fluxpool.hpp:

inline void fluxpool::set(double v, unit_types u, bool track = false,
                          string pool_name = "test") {
  name = pool_name;
  if (v < 0) {
    H_ASSERT(v >= 0, "Flux and pool values may not be negative in 2 " + name);
  }

even updating the default pool string name from “?”, rebuilding hector, doing devtools::load_all here so it’s updated, the resulting error is still “Flux and pool values may not be negative in 2 ?”
added the 2 to the error statement so I knew Rstudio was seeing update in C++ code’.

Solution Part 1: only use an NBP constraint with positive values.

In practice, this just means only using years >1970 or so as the constraint.
for GCP, depending on how you smooth it, this can be 1972. But if you try to use 1972 on with hector data, need to figure out a different diffusivity. very dumb.
This led to some different errors, sometimes, that made me think I was on the right track but still ultimately gave the same flux values may not be negative error most of the time.

Solution Part 2: Lower the ocean heat diffusivity from its default value of 2.38 cm2/s to 2.31 cm2/s (for a 245).

ACS had run some NBP constraints a year+ ago and they were NOT as finicky to get to solve, but that version of hector had different beta, q10, and diffusivity.
the old (1.16) vs new value of diffusivity was the biggest difference (new was about double old), so started at 2.38 and just started incremementing lower until the constraint solved.
definitely had some runs with 2.34 that solved but can’t replicate. I think this is not a real instability in hector code but is instead user error - in rstudio, you can’t cmd+enter a bunch of rhector function calls together - some end up kind of skipped because the newcore and setvar functions seem to take an extra moment and need to be run on their own before proceeding?
but also only sometimes. Seems to always be fine when knitting at least.

functions

handful of steps with the above, so make a function instead of copying and pasting every time.
adjusting ocean heat diffusivity to 2.31 cm2/s

adjust_hector_run <- function(initial_ini,
                              scen_name,
                              diffusivity = 2.31){
    
    rcp <- initial_ini
    scenario_file <- paste0("input/hector_",rcp,".ini")
    ini_file <- system.file(scenario_file, package="hector") 
    
    core <- hector::newcore(ini_file, suppresslogging = FALSE, name=scen_name)
    hector::run(core,runtodate = 2100)
    
    setvar(core, NA, DIFFUSIVITY(), diffusivity, "cm2/s") # default 2.38
    print(fetchvars(core, NA, DIFFUSIVITY()))
    
    hector::run(core,runtodate = 2100)
    
    return(core)
}

output_extractor <- function(core){
    fetchvars(core,
              dates=1750:2100,
              vars=c(GLOBAL_TAS(),
                     CONCENTRATIONS_CO2(), 
                     NBP(),
                     LL_OCEAN_UPTAKE(),
                     HL_OCEAN_UPTAKE(),
                     SOIL_C(),
                     VEG_C()))  %>%
        mutate(variable = if_else( variable == 'NBP', 'A.NBP',
                               if_else(variable == 'CO2_concentration', 'B.CO2_concentration',
                                       if_else(variable == 'global_tas', 'C.global_tas',
                                               if_else(variable == 'LL_ocean_uptake', 'F.LL_ocean_uptake',
                                                       if_else(variable == 'HL_ocean_uptake', 'G.HL_ocean_uptake',
                                                               if_else( variable == 'soil_c', 'D.soil_c',
                                                                        if_else(variable == 'veg_c', 'E.veg_c',
                                                                                variable)))))))) %>%
    mutate(variable = paste0(variable, '~', units)) 
    
}

Check effect of ocean param change

# Default hector
rcp <- "ssp245"
scenario_file <- paste0("input/hector_",rcp,".ini")
ini_file245 <- system.file(scenario_file, package="hector") 

core245 <- hector::newcore(ini_file245, name = 'ssp245_default')
hector::run(core245,runtodate = 2100)

## Hector core: ssp245_default
## Start date:  1745
## End date:    2300
## Current date:    2100
## Input file:  /Users/snyd535/Documents/GitHub/hector/inst/input/hector_ssp245.ini

# Core with adjusted params that work for constraint
#reset(core_test)
core_test <- adjust_hector_run(initial_ini = 'ssp245',
                                scen_name = 'ssp245_diff2.31')

##          scenario year variable value units
## 1 ssp245_diff2.31   NA     diff  2.31 cm2/s

## Auto-resetting core to 0

hector::run(core_test,runtodate = 2100)

## Hector core: ssp245_diff2.31
## Start date:  1745
## End date:    2300
## Current date:    2100
## Input file:  /Users/snyd535/Documents/GitHub/hector/inst/input/hector_ssp245.ini

# extract outputs
bind_rows(output_extractor(core245),
          output_extractor(core_test)) ->
    plot_out

ggplot(plot_out) +
    geom_line(mapping = aes(x=year, y = value, color = scenario, linetype = scenario)) +
    facet_wrap(~variable, scales = 'free_y', nrow=2) +
    ggtitle('Effect of diffusivity 2.31 cm2/s (old = 2.38)')

note visible HL_ocean_uptake potential instability in ~2070, otherwise differences are:

plot_out %>%
    filter(scenario == 'ssp245_default') %>%
    rename(ssp245_default = value) %>%
    select(-scenario) %>%
    left_join(plot_out %>%
                  filter(scenario == 'ssp245_diff2.31') %>%
                  rename(ssp245_diff2.31 = value) %>%
                  select(-scenario),
              by = c('year', 'variable', 'units') ) %>%
    group_by(variable, units) %>%
    summarize(maxdiff = max(abs(ssp245_default-ssp245_diff2.31)),
              maxnormdiff = max(abs(ssp245_default-ssp245_diff2.31)) /max(abs(ssp245_default))) %>%
    ungroup %>%
    knitr::kable()

## `summarise()` has grouped output by 'variable'. You can override using the
## `.groups` argument.

variable	units	maxdiff	maxnormdiff
A.NBP~Pg C/yr	Pg C/yr	0.0076322	0.0024744
B.CO2_concentration~ppmv CO2	ppmv CO2	0.3029845	0.0005408
C.global_tas~degC	degC	0.0132185	0.0048332
D.soil_c~Pg C	Pg C	0.4119093	0.0001965
E.veg_c~Pg C	Pg C	0.1063067	0.0001582
F.LL_ocean_uptake~Pg C	Pg C	0.0025119	0.0009716
G.HL_ocean_uptake~Pg C	Pg C	0.0005476	0.0002809

take a look at ocean in more detail

ggplot(plot_out %>% filter(grepl('ocean', variable))) +
    geom_line(mapping = aes(x=year, y = value, color = scenario, linetype = scenario)) +
    facet_wrap(~variable, scales = 'free_y', nrow=2) +
    ggtitle('Effect of diffusivity 2.31 cm2/s (old = 2.38)')

So the potential instability in the ocean uptake is present for both values of diffusivity

So not surprisingly based on a lot of other things we have seen over the last couple years with hector, the ocean is the part that needs maintenance.

good to know diffusivity is a lever. Will do some sensitivity tests further down sweeping out different diffusivity values vs different magnitudes of NBP (will explain more there).
zoom in some more on terrestrial carbon quantities

ggplot(plot_out %>% filter(grepl('NBP', variable) | grepl('veg', variable) | grepl('soil', variable)) )+
    geom_line(mapping = aes(x=year, y = value, color = scenario, linetype = scenario)) +
    facet_wrap(~variable, scales = 'free_y', nrow=3) +
    ggtitle('Effect of diffusivity 2.31 cm2/s (old = 2.38)')

nothing obvious stands out.

So from this point, all Hector runs will be done with 2.31cm2/s for diffusivity so the NBP constraints will actually run without running into negative flux pools.

0-level test: Hector 245 1970-2100 constraint

If we extract the entire nbp trajectory from a hector run and make it the constraint for a run with otherwise identical settings, do we get identical results?

# clean up some variables
rm(core_test)
rm(core245)
rm(plot_out)

# initialize a holder of outputs 
output_holder <- list()

# Do the run to generate the constraing
core <- adjust_hector_run(initial_ini = 'ssp245',
                                scen_name = 'ssp245_gen')

##     scenario year variable value units
## 1 ssp245_gen   NA     diff  2.31 cm2/s

## Auto-resetting core to 0

output_holder[['ssp245_gen']] <- output_extractor(core)
rm(core)


# Extract Constraint 
output_holder[['ssp245_gen']] %>% 
    filter(grepl('NBP', variable),
           year >= 1970,
           year<= 2100) ->
    nbp_constrainer

# Set up core for constraint, set constraint, run wtih constraint
core <- adjust_hector_run(initial_ini = 'ssp245',
                          scen_name = 'ssp245_1970~2100')

##           scenario year variable value units
## 1 ssp245_1970~2100   NA     diff  2.31 cm2/s

## Auto-resetting core to 0

setvar(core, nbp_constrainer$year, NBP_CONSTRAIN(), nbp_constrainer$value, "Pg C/yr")

hector::run(core,runtodate = 2100)

## Auto-resetting core to 1969

## Hector core: ssp245_1970~2100
## Start date:  1745
## End date:    2300
## Current date:    2100
## Input file:  /Users/snyd535/Documents/GitHub/hector/inst/input/hector_ssp245.ini

output_holder[['ssp245_1970~2100']] <- output_extractor(core)
rm(core)

ggplot(do.call(rbind, output_holder) %>% filter(grepl('245', scenario))) +
    geom_line(mapping = aes(x=year, y = value, color = scenario, linetype = scenario)) +
    facet_wrap(~variable, scales = 'free_y', nrow=2) +
    ggtitle('testing NBP constraint SSP245')

Previously only compared NBP before and after applying constraint.
clear from looking at these additional quantities, it’s not working correctly.
differences in soil and veg carbon are possibly a bit weird (depending on how constraint implemented in code). But in theory, calculations from the land box don’t go anywhere when there is an NBP constraint, so a difference is ok as long as the land calculations don’t co anywhere.
But the ocean uptakes differing as well suggests the land calculations ARE going somewhere and being used instead of the constraint.
HYPOTHESIS in code, the NBP constraint is overwriting some values in some places (which is why NBP constraint returns identical) but probably not in all of the places it should.
So constraint being on = land box evolves differently, with lower soil and veg carbon stocks. Lower land carbon stocks (soil especially) means lower ocean uptake in reaction? is that right directions?

## `summarise()` has grouped output by 'variable'. You can override using the
## `.groups` argument.

SSP245 NBP constraint testing
variable	units	maxdiff	maxnormdiff
A.NBP~Pg C/yr	Pg C/yr	0.0000000	0.0000000
B.CO2_concentration~ppmv CO2	ppmv CO2	4.2906869	0.0076546
C.global_tas~degC	degC	0.0209422	0.0076212
D.soil_c~Pg C	Pg C	52.2774744	0.0249493
E.veg_c~Pg C	Pg C	8.2995916	0.0123482
F.LL_ocean_uptake~Pg C	Pg C	0.1860536	0.0719548
G.HL_ocean_uptake~Pg C	Pg C	0.1499027	0.0768922

0-level test: Hector 126 1970-2100 constraint

Know ssp245 can be weird (github issue), so double check other scenarios

# Do the run to generate the constraing
core <- adjust_hector_run(initial_ini = 'ssp126',
                                scen_name = 'ssp126_gen')

##     scenario year variable value units
## 1 ssp126_gen   NA     diff  2.31 cm2/s

## Auto-resetting core to 0

output_holder[['ssp126_gen']] <- output_extractor(core)
rm(core)


# Extract Constraint 
output_holder[['ssp126_gen']] %>% 
    filter(grepl('NBP', variable),
           year >= 1970,
           year<= 2100) ->
    nbp_constrainer

# Set up core for constraint, set constraint, run wtih constraint
core <- adjust_hector_run(initial_ini = 'ssp126',
                          scen_name = 'ssp126_1970~2100')

##           scenario year variable value units
## 1 ssp126_1970~2100   NA     diff  2.31 cm2/s

## Auto-resetting core to 0

setvar(core, nbp_constrainer$year, NBP_CONSTRAIN(), nbp_constrainer$value, "Pg C/yr")

hector::run(core,runtodate = 2100)

## Auto-resetting core to 1969

## Hector core: ssp126_1970~2100
## Start date:  1745
## End date:    2300
## Current date:    2100
## Input file:  /Users/snyd535/Documents/GitHub/hector/inst/input/hector_ssp126.ini

output_holder[['ssp126_1970~2100']] <- output_extractor(core)
rm(core)


ggplot(do.call(rbind, output_holder) %>% filter(grepl('126', scenario))) +
    geom_line(mapping = aes(x=year, y = value, color = scenario, linetype = scenario)) +
    facet_wrap(~variable, scales = 'free_y', nrow=2) +
    ggtitle('testing NBP constraint SSP126')

so not just a 245 quirk

0-level test: Hector 585 1970-2100 constraint

# Do the run to generate the constraing
core <- adjust_hector_run(initial_ini = 'ssp585',
                                scen_name = 'ssp585_gen')

##     scenario year variable value units
## 1 ssp585_gen   NA     diff  2.31 cm2/s

## Auto-resetting core to 0

output_holder[['ssp585_gen']] <- output_extractor(core)
rm(core)


# Extract Constraint 
output_holder[['ssp585_gen']] %>% 
    filter(grepl('NBP', variable),
           year >= 1970,
           year<= 2100) ->
    nbp_constrainer

# Set up core for constraint, set constraint, run wtih constraint
core <- adjust_hector_run(initial_ini = 'ssp585',
                          scen_name = 'ssp585_1970~2100')

##           scenario year variable value units
## 1 ssp585_1970~2100   NA     diff  2.31 cm2/s

## Auto-resetting core to 0

setvar(core, nbp_constrainer$year, NBP_CONSTRAIN(), nbp_constrainer$value, "Pg C/yr")

hector::run(core,runtodate = 2100)

## Auto-resetting core to 1969

## Hector core: ssp585_1970~2100
## Start date:  1745
## End date:    2300
## Current date:    2100
## Input file:  /Users/snyd535/Documents/GitHub/hector/inst/input/hector_ssp585.ini

output_holder[['ssp585_1970~2100']] <- output_extractor(core)
rm(core)


ggplot(do.call(rbind, output_holder) %>% filter(grepl('585', scenario))) +
    geom_line(mapping = aes(x=year, y = value, color = scenario, linetype = scenario)) +
    facet_wrap(~variable, scales = 'free_y', nrow=2) +
    ggtitle('testing NBP constraint SSP585')

More constraints

hector 245 1970-2021 constraint

Seems useful to understand what happens when you just constrain the (usable) historical period.
usable in the sense that need NBP constraint values not to be negative, so years > 1970.
cutting at 2021 because that’s the last year in the GCP data I have.

# Constraint with updated params
output_holder[['ssp245_gen']] %>% 
    filter(grepl('NBP', variable),
           year >= 1970,
           year<= 2021) ->
    nbp_constrainer

# Set up core for constraint, set constraint, run wtih constraint
core <- adjust_hector_run(initial_ini = 'ssp245',
                          scen_name = 'ssp245_1970~2021')

##           scenario year variable value units
## 1 ssp245_1970~2021   NA     diff  2.31 cm2/s

## Auto-resetting core to 0

setvar(core, nbp_constrainer$year, NBP_CONSTRAIN(), nbp_constrainer$value, "Pg C/yr")

hector::run(core,runtodate = 2100)

## Auto-resetting core to 1969

## Hector core: ssp245_1970~2021
## Start date:  1745
## End date:    2300
## Current date:    2100
## Input file:  /Users/snyd535/Documents/GitHub/hector/inst/input/hector_ssp245.ini

output_holder[['ssp245_1970~2021']] <- output_extractor(core)
rm(core)


ggplot(do.call(rbind, output_holder) %>% filter(grepl('245', scenario))) +
    geom_line(mapping = aes(x=year, y = value, color = scenario, linetype = scenario)) +
    facet_wrap(~variable, scales = 'free_y', nrow=2) +
    geom_vline(xintercept = 2021, size = 0.2)+
    ggtitle('testing NBP constraint SSP245')

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

By terminating the constraint in 2021, all variables are able to drift closer to their default hector behavior.
NOTE the abrupt jump in NBP when the constraint ceases. HYPOTHESIS The background land calculations going on while the constraint is in place are immediately what kicks in at the end of the constraint, rather than the constraint value.
In 2020, background calculation produces lower carbon stocks in land = abrupt jump in land sink? can’t tell if that’s right but blue vs green and red vs green curve comparison in this figure is the key to figuring it out, I think.

GCP NBP constraint 1972-2021

A few years in the GCP data from 1970-1980 have nbp<0.

This seems to cause the NBP constraint Hector actually sees to remain at 0, leading to a huge rebound spike when the constraint ends in 2021.

Just droping the specific years with negative nbp causes Hector to fill in 0s for those years, leading to a rebound at each of these years (not shown).

# Run with GCP NBP
gcp_data %>% 
    filter(year >= 1972,
           year<= 2021) %>%
    #filter(nbp>0) %>%
    mutate(value = 0.001*nbp,
           Units = "Pg C/yr")->
    nbp_constrainer

# Set up core for constraint, set constraint, run wtih constraint
core <- adjust_hector_run(initial_ini = 'ssp245',
                          scen_name = 'ssp245_GCP1972~2021')

##              scenario year variable value units
## 1 ssp245_GCP1972~2021   NA     diff  2.31 cm2/s

## Auto-resetting core to 0

setvar(core, nbp_constrainer$year, NBP_CONSTRAIN(), nbp_constrainer$value, "Pg C/yr")


hector::run(core,runtodate = 2100)

## Auto-resetting core to 1971

## Hector core: ssp245_GCP1972~2021
## Start date:  1745
## End date:    2300
## Current date:    2100
## Input file:  /Users/snyd535/Documents/GitHub/hector/inst/input/hector_ssp245.ini

output_holder[['ssp245_GCP1972~2021']] <- output_extractor(core)
rm(core)


ggplot(do.call(rbind, output_holder) %>% filter(grepl('245', scenario))) +
    geom_line(mapping = aes(x=year, y = value, color = scenario, linetype = scenario)) +
    facet_wrap(~variable, scales = 'free_y', nrow=2) +
    geom_vline(xintercept = 2021, size = 0.2)+
    ggtitle('testing NBP constraint SSP245')

Near end of GCP constaint (~2015), NBP is almost twice as large as in the hector default. This leads to a rebound dip.

245 runs magnified

ggplot(do.call(rbind, output_holder) %>% filter(grepl('245', scenario), year >= 1965, year <=2050)) +
    geom_line(mapping = aes(x=year, y = value, color = scenario, linetype = scenario)) +
    facet_wrap(~variable, scales = 'free_y', nrow=2) +
    geom_vline(xintercept = 2021, size = 0.2)+ 
    ggtitle('testing NBP constraint SSP245')

hector 585 1970-2021 constraint

# Constraint with updated params
output_holder[['ssp585_gen']] %>% 
    filter(grepl('NBP', variable),
           year >= 1970,
           year<= 2021) ->
    nbp_constrainer

# Set up core for constraint, set constraint, run wtih constraint
core <- adjust_hector_run(initial_ini = 'ssp585',
                          scen_name = 'ssp585_1970~2021')

##           scenario year variable value units
## 1 ssp585_1970~2021   NA     diff  2.31 cm2/s

## Auto-resetting core to 0

setvar(core, nbp_constrainer$year, NBP_CONSTRAIN(), nbp_constrainer$value, "Pg C/yr")

hector::run(core,runtodate = 2100)

## Auto-resetting core to 1969

## Hector core: ssp585_1970~2021
## Start date:  1745
## End date:    2300
## Current date:    2100
## Input file:  /Users/snyd535/Documents/GitHub/hector/inst/input/hector_ssp585.ini

output_holder[['ssp585_1970~2021']] <- output_extractor(core)
rm(core)


ggplot(do.call(rbind, output_holder) %>% filter(grepl('585', scenario))) +
    geom_line(mapping = aes(x=year, y = value, color = scenario, linetype = scenario)) +
    facet_wrap(~variable, scales = 'free_y', nrow=2) +
    geom_vline(xintercept = 2021, size = 0.2)+
    ggtitle('testing NBP constraint SSP585')

GCP 585 1972-2021

# Run with GCP NBP
gcp_data %>% 
    filter(year >= 1972,
           year<= 2021) %>%
    #filter(nbp>0) %>%
    mutate(value = 0.001*nbp,
           Units = "Pg C/yr")->
    nbp_constrainer

# Set up core for constraint, set constraint, run wtih constraint
core <- adjust_hector_run(initial_ini = 'ssp585',
                          scen_name = 'ssp585_GCP1972~2021')

##              scenario year variable value units
## 1 ssp585_GCP1972~2021   NA     diff  2.31 cm2/s

## Auto-resetting core to 0

setvar(core, nbp_constrainer$year, NBP_CONSTRAIN(), nbp_constrainer$value, "Pg C/yr")


hector::run(core,runtodate = 2100)

## Auto-resetting core to 1971

## Hector core: ssp585_GCP1972~2021
## Start date:  1745
## End date:    2300
## Current date:    2100
## Input file:  /Users/snyd535/Documents/GitHub/hector/inst/input/hector_ssp585.ini

output_holder[['ssp585_GCP1972~2021']] <- output_extractor(core)
rm(core)


ggplot(do.call(rbind, output_holder) %>% filter(grepl('585', scenario))) +
    geom_line(mapping = aes(x=year, y = value, color = scenario, linetype = scenario)) +
    facet_wrap(~variable, scales = 'free_y', nrow=2) +
    geom_vline(xintercept = 2021, size = 0.2)+
    ggtitle('testing NBP constraint SSP585')

GCP 585 1972-2000

# Run with GCP NBP
gcp_data %>% 
    filter(year >= 1972,
           year<= 2000) %>%
    #filter(nbp>0) %>%
    mutate(value = 0.001*nbp,
           Units = "Pg C/yr")->
    nbp_constrainer

# Set up core for constraint, set constraint, run wtih constraint
core <- adjust_hector_run(initial_ini = 'ssp585',
                          scen_name = 'ssp585_GCP1972~2000')

##              scenario year variable value units
## 1 ssp585_GCP1972~2000   NA     diff  2.31 cm2/s

## Auto-resetting core to 0

setvar(core, nbp_constrainer$year, NBP_CONSTRAIN(), nbp_constrainer$value, "Pg C/yr")


hector::run(core,runtodate = 2100)

## Auto-resetting core to 1971

## Hector core: ssp585_GCP1972~2000
## Start date:  1745
## End date:    2300
## Current date:    2100
## Input file:  /Users/snyd535/Documents/GitHub/hector/inst/input/hector_ssp585.ini

output_holder[['ssp585_GCP1972~2000']] <- output_extractor(core)
rm(core)


ggplot(do.call(rbind, output_holder) %>% filter(grepl('585', scenario))) +
    geom_line(mapping = aes(x=year, y = value, color = scenario, linetype = scenario)) +
    facet_wrap(~variable, scales = 'free_y', nrow=2) +
    geom_vline(xintercept = 2021, size = 0.2)+ 
    geom_vline(xintercept = 2021, size = 0.1, color = 'grey')+ 
    ggtitle('testing NBP constraint SSP585')

585 runs magnified

ggplot(do.call(rbind, output_holder) %>% filter(grepl('585', scenario), year >= 1965, year <=2050)) +
    geom_line(mapping = aes(x=year, y = value, color = scenario, linetype = scenario)) +
    facet_wrap(~variable, scales = 'free_y', nrow=2) +
    geom_vline(xintercept = 2021, size = 0.2)+ 
    geom_vline(xintercept = 2021, size = 0.1, color = 'grey')+ 
    ggtitle('testing NBP constraint SSP585')

126 1970-2021

Constraints across scenarios

# Constraint with updated params
output_holder[['ssp126_gen']] %>% 
    filter(grepl('NBP', variable),
           year >= 1970,
           year<= 2021) ->
    nbp_constrainer

# Set up core for constraint, set constraint, run wtih constraint
core <- adjust_hector_run(initial_ini = 'ssp126',
                          scen_name = 'ssp126_1970~2021')

##           scenario year variable value units
## 1 ssp126_1970~2021   NA     diff  2.31 cm2/s

## Auto-resetting core to 0

setvar(core, nbp_constrainer$year, NBP_CONSTRAIN(), nbp_constrainer$value, "Pg C/yr")

hector::run(core,runtodate = 2100)

## Auto-resetting core to 1969

## Hector core: ssp126_1970~2021
## Start date:  1745
## End date:    2300
## Current date:    2100
## Input file:  /Users/snyd535/Documents/GitHub/hector/inst/input/hector_ssp126.ini

output_holder[['ssp126_1970~2021']] <- output_extractor(core)
rm(core)


ggplot(do.call(rbind, output_holder) %>% filter(!grepl('GCP', scenario), year >= 1965, year <=2050)) +
    geom_line(mapping = aes(x=year, y = value, color = scenario, linetype = scenario)) +
    facet_wrap(~variable, scales = 'free_y', nrow=2) +
    geom_vline(xintercept = 2021, size = 0.2)+
    ggtitle('testing NBP constraint')

Conclusions

NBP constraint is not properly passing/talking to all other parts of Hector.

So my guess of what’s happening, just based on the plots, is the land calculations still run in the background and the NBP constraint is supposed to just ignore them. But somewhere, the land calculations aren’t actually being ignored and the ocean uptake reacts to them.
So I think Atm C is getting the nbp constraint for land, but also the ocean differences reacting to background land caluclations insteadof NBP and that’s what makes AtmC a little different?
No idea:

A. Where in the code the land calculations that should be ignored are getting passed on to anything

B. Why those background land calculations are different

There is also:

HL ocean uptake plots also suggest potential for numerical instability going on there.

DISREGARD from here other tests

Can update the below when NBP constraint actually mechanically works.

all of above was just trying to get NBP cosntraint to run at all.

Need:

to make sure that for an ordered set of NBP constraints, T and CO2 should be following sensible order
sensitivity on magnitude of constraint vs diffusivity

Ordered constraints 1970-2021 and T, CO2

keep it simple - couple multiples of 245 NBP with 2.31

# Get constraint:
core_constrain <- hector::newcore(ini_file245, name = '245_diff2.31_generate')
setvar(core_constrain, NA, DIFFUSIVITY(), 2.31, "cm2/s") # default 2.38


hector::run(core_constrain,runtodate = 2050)

nbp <- fetchvars(core_constrain, 1965:2050, NBP())

out  <-fetchvars(core_constrain,dates=1965:2050,vars=c(GLOBAL_TAS(),
                                                       CONCENTRATIONS_CO2(), 
                                                       NBP(),
                                                       LL_OCEAN_UPTAKE(),
                                                       HL_OCEAN_UPTAKE())) %>%
    mutate(k=1)


# Constraint with updated params
nbp %>% 
    filter(scenario == '245_diff2.31_generate',
           year >= 1970,
           year<= 2021) ->
    nbp_constrainer


for (k in c(0.5, 1, 2, 5, 8)){
    core_k <- hector::newcore(ini_file245, name = paste0('245_diff2.31_constrain_k',k))
    setvar(core_k, NA, DIFFUSIVITY(), 2.31, "cm2/s")
    
    setvar(core_k, nbp_constrainer$year, NBP_CONSTRAIN(), k*nbp_constrainer$value, "Pg C/yr")
    
    hector::run(core_k,runtodate = 2050)

    
    fetchvars(core_k,dates=1965:2050,vars=c(GLOBAL_TAS(),
                                            CONCENTRATIONS_CO2(), 
                                            NBP(),
                                            LL_OCEAN_UPTAKE(),
                                            HL_OCEAN_UPTAKE())) %>%
        mutate(k=k) %>%
        bind_rows(out) -> 
        out
    rm(core_k)
}


# Order the variables for faceting
out %>%
    mutate(variable = if_else( variable == 'NBP', 'A.NBP',
                               if_else(variable == 'CO2_concentration', 'B.CO2_concentration',
                                       if_else(variable == 'global_tas', 'C.global_tas',
                                               if_else(variable == 'LL_ocean_uptake', 'D.LL_ocean_uptake',
                                                       if_else(variable == 'HL_ocean_uptake', 'E.HL_ocean_uptake',
                                                               variable)))))) %>%
    mutate(variable = paste0(variable, '~', units)) ->
    out2

ggplot(out2) +
    geom_line(mapping = aes(x=year, y = value, color = k, linetype = scenario)) +
    facet_wrap(~variable, scales = 'free_y', nrow=1) +
    geom_vline(xintercept = 2021, size = 0.3) +
    scale_color_continuous() +
    ggtitle('diffusivity 2.31 cm2/s')

Ordered constraints 1970-2050 and T, CO2

keep it simple - couple multiples of 245 NBP with 2.31

# Get constraint:
nbp <- fetchvars(core_constrain, 1965:2050, NBP())

out  <-fetchvars(core_constrain,dates=1965:2050,vars=c(GLOBAL_TAS(),
                                                       CONCENTRATIONS_CO2(), 
                                                       NBP(),
                                                       LL_OCEAN_UPTAKE(),
                                                       HL_OCEAN_UPTAKE())) %>%
    mutate(k=1)


# Constraint with updated params
nbp %>% 
    filter(scenario == '245_diff2.31_generate',
           year >= 1970,
           year<= 2050) ->
    nbp_constrainer


for (k in c(0.5, 1, 2, 5)){
    core_k <- hector::newcore(ini_file245, name = paste0('245_diff2.31_constrain_k',k))
    setvar(core_k, NA, DIFFUSIVITY(), 2.31, "cm2/s")
    
    setvar(core_k, nbp_constrainer$year, NBP_CONSTRAIN(), k*nbp_constrainer$value, "Pg C/yr")
    
    hector::run(core_k,runtodate = 2050)

    
    fetchvars(core_k,dates=1965:2050,vars=c(GLOBAL_TAS(),
                                            CONCENTRATIONS_CO2(), 
                                            NBP(),
                                            LL_OCEAN_UPTAKE(),
                                            HL_OCEAN_UPTAKE())) %>%
        mutate(k=k) %>%
        bind_rows(out) -> 
        out
    rm(core_k)
}


# Order the variables for faceting
out %>%
    mutate(variable = if_else( variable == 'NBP', 'A.NBP',
                               if_else(variable == 'CO2_concentration', 'B.CO2_concentration',
                                       if_else(variable == 'global_tas', 'C.global_tas',
                                               if_else(variable == 'LL_ocean_uptake', 'D.LL_ocean_uptake',
                                                       if_else(variable == 'HL_ocean_uptake', 'E.HL_ocean_uptake',
                                                               variable)))))) %>%
    mutate(variable = paste0(variable, '~', units)) ->
    out2

ggplot(out2) +
    geom_line(mapping = aes(x=year, y = value, color = k, linetype = scenario)) +
    facet_wrap(~variable, scales = 'free_y', nrow=1) +
    # geom_vline(xintercept = 2021, size = 0.3) +
    scale_color_continuous() +
    ggtitle('diffusivity 2.31 cm2/s')

lower diffusivity 1970-2021

use 1.16, show that nbp/co2/temp don’t behave the way they should.

# Get constraint:
core_constrain <- hector::newcore(ini_file245, name = '245_diff1.16_generate')
setvar(core_constrain, NA, DIFFUSIVITY(), 1.16, "cm2/s") # default 2.38


hector::run(core_constrain,runtodate = 2050)

nbp <- fetchvars(core_constrain, 1965:2050, NBP())

out  <-fetchvars(core_constrain,dates=1965:2050,vars=c(GLOBAL_TAS(),
                                                       CONCENTRATIONS_CO2(), 
                                                       NBP(),
                                                       LL_OCEAN_UPTAKE(),
                                                       HL_OCEAN_UPTAKE())) %>%
    mutate(k=1)


# Constraint with updated params
nbp %>% 
    filter(scenario == '245_diff1.16_generate',
           year >= 1970,
           year<= 2021) ->
    nbp_constrainer


for (k in c(0.5, 1, 2, 5, 8)){
    core_k <- hector::newcore(ini_file245, name = paste0('245_diff1.16_constrain_k',k))
    setvar(core_k, NA, DIFFUSIVITY(), 1.16, "cm2/s")
    
    setvar(core_k, nbp_constrainer$year, NBP_CONSTRAIN(), k*nbp_constrainer$value, "Pg C/yr")
    
    hector::run(core_k,runtodate = 2050)

    
    fetchvars(core_k,dates=1965:2050,vars=c(GLOBAL_TAS(),
                                            CONCENTRATIONS_CO2(), 
                                            NBP(),
                                            LL_OCEAN_UPTAKE(),
                                            HL_OCEAN_UPTAKE())) %>%
        mutate(k=k) %>%
        bind_rows(out) -> 
        out
    rm(core_k)
}


# Order the variables for faceting
out %>%
    mutate(variable = if_else( variable == 'NBP', 'A.NBP',
                               if_else(variable == 'CO2_concentration', 'B.CO2_concentration',
                                       if_else(variable == 'global_tas', 'C.global_tas',
                                               if_else(variable == 'LL_ocean_uptake', 'D.LL_ocean_uptake',
                                                       if_else(variable == 'HL_ocean_uptake', 'E.HL_ocean_uptake',
                                                               variable)))))) %>%
    mutate(variable = paste0(variable, '~', units)) ->
    out2

ggplot(out2) +
    geom_line(mapping = aes(x=year, y = value, color = k, linetype = scenario)) +
    facet_wrap(~variable, scales = 'free_y', nrow=1) +
    geom_vline(xintercept = 2021, size = 0.3) +
    scale_color_continuous() +
    ggtitle('diffusivity 1.16 cm2/s')

lower diffusivity 1970-2050

use 1.16, show that nbp/co2/temp don’t behave the way they should.

# Get constraint:
nbp <- fetchvars(core_constrain, 1965:2050, NBP())

out  <-fetchvars(core_constrain,dates=1965:2050,vars=c(GLOBAL_TAS(),
                                                       CONCENTRATIONS_CO2(), 
                                                       NBP(),
                                                       LL_OCEAN_UPTAKE(),
                                                       HL_OCEAN_UPTAKE())) %>%
    mutate(k=1)


# Constraint with updated params
nbp %>% 
    filter(scenario == '245_diff1.16_generate',
           year >= 1970,
           year<= 2050) ->
    nbp_constrainer


for (k in c(0.5, 1, 2, 5)){
    core_k <- hector::newcore(ini_file245, name = paste0('245_diff1.16_constrain_k',k))
    setvar(core_k, NA, DIFFUSIVITY(), 1.16, "cm2/s")
    
    setvar(core_k, nbp_constrainer$year, NBP_CONSTRAIN(), k*nbp_constrainer$value, "Pg C/yr")
    
    hector::run(core_k,runtodate = 2050)

    
    fetchvars(core_k,dates=1965:2050,vars=c(GLOBAL_TAS(),
                                            CONCENTRATIONS_CO2(), 
                                            NBP(),
                                            LL_OCEAN_UPTAKE(),
                                            HL_OCEAN_UPTAKE())) %>%
        mutate(k=k) %>%
        bind_rows(out) -> 
        out
    rm(core_k)
}


# Order the variables for faceting
out %>%
    mutate(variable = if_else( variable == 'NBP', 'A.NBP',
                               if_else(variable == 'CO2_concentration', 'B.CO2_concentration',
                                       if_else(variable == 'global_tas', 'C.global_tas',
                                               if_else(variable == 'LL_ocean_uptake', 'D.LL_ocean_uptake',
                                                       if_else(variable == 'HL_ocean_uptake', 'E.HL_ocean_uptake',
                                                               variable)))))) %>%
    mutate(variable = paste0(variable, '~', units)) ->
    out2

ggplot(out2) +
    geom_line(mapping = aes(x=year, y = value, color = k, linetype = scenario)) +
    facet_wrap(~variable, scales = 'free_y', nrow=1) +
    geom_vline(xintercept = 2021, size = 0.3) +
    scale_color_continuous() +
    ggtitle('diffusivity 1.16 cm2/s')

replicate NBP-T-C issue

# GCP
gcp_data %>% 
    filter(year >= 1972,
           year<= 2021) %>%
    #filter(nbp>0) %>%
    mutate(nbp = 0.001*nbp,
           Units = "Pg C/yr")->
    nbp_constrainer


coreGCP <- hector::newcore(ini_file245, name = 'GCP_diff2.31_constraint_1972~2021')
# setvar(coreGCP, NA, BETA(), 0.55, "(unitless)")
# setvar(coreGCP, NA, Q10_RH(), 2.2, "(unitless)") 
setvar(coreGCP, NA, DIFFUSIVITY(), 1.16, "cm2/s") 

setvar(coreGCP,nbp_constrainer$year,"NBP_constrain",nbp_constrainer$nbp,"Pg C/yr")

hector::run(coreGCP,runtodate = 2050)


out  <-fetchvars(coreGCP, dates=1965:2050,vars=c(GLOBAL_TAS(),
                                                       CONCENTRATIONS_CO2(), 
                                                       NBP(),
                                                       LL_OCEAN_UPTAKE(),
                                                       HL_OCEAN_UPTAKE())) %>%
    mutate(k=1)


for (k in c(2, 6)){
    core_k <- hector::newcore(ini_file245, name = paste0('245_k',k))
    # setvar(coreGCP, NA, BETA(), 0.55, "(unitless)") # have to be consistent with above
    # setvar(coreGCP, NA, Q10_RH(),2.2, "(unitless)") # have to be consistent with above
    setvar(core_k, NA, DIFFUSIVITY(), 1.16, "cm2/s")
    
    setvar(core_k, nbp_constrainer$year, NBP_CONSTRAIN(), k*nbp_constrainer$nbp, "Pg C/yr")
    
    hector::run(core_k,runtodate = 2050)

    
    fetchvars(core_k,dates=1965:2050,vars=c(GLOBAL_TAS(),
                                            CONCENTRATIONS_CO2(), 
                                            NBP(),
                                            LL_OCEAN_UPTAKE(),
                                            HL_OCEAN_UPTAKE())) %>%
        mutate(k=k) %>%
        bind_rows(out) -> 
        out
    rm(core_k)
}


# Order the variables for faceting
out %>%
    mutate(variable = if_else( variable == 'NBP', 'A.NBP',
                               if_else(variable == 'CO2_concentration', 'B.CO2_concentration',
                                       if_else(variable == 'global_tas', 'C.global_tas',
                                               if_else(variable == 'LL_ocean_uptake', 'D.LL_ocean_uptake',
                                                       if_else(variable == 'HL_ocean_uptake', 'E.HL_ocean_uptake',
                                                               variable)))))) %>%
    mutate(variable = paste0(variable, '~', units)) ->
    out2

ggplot(out2 %>% filter(year <= 2000)) +
    geom_line(mapping = aes(x=year, y = value, color = k, linetype = scenario)) +
    facet_wrap(~variable, scales = 'free_y', nrow=1) +
    geom_vline(xintercept = 1994, size = 0.3) +
    scale_color_continuous() #+
    #ggtitle('diffusivity 1.16 cm2/s, beta=0.55, q10=2.2 \n just changing diff is not enough')

diff: 1.16, beta:.55, q10: 2.2

can’t replicate but ocean uptake in 1993-1995 interesting: in default, both LL and HL ocean are a valley at the same time as NBP valley.

But when that NBP valley is made more extreme, HL and LL become peaks

diff:2.31, beta:0.53, q10:1.76

TODO have plots for both parameterization, not just toggle in code

ALSO TODO: why changing beta and q10 changing plots for same diffusivity (double check side by side) - have to use diff 1.16 to solve with both sets of beta/q10 but identical as should be (think)

make a copy of the old notebook and see if can re-generate there with old params and new hector.

sweep diffusivity vs gcp magnitude

test doing one year at a time (like DW work)

is there a rebound when feeding on next time step

ALSO double check gcamlandc - is it passing emissions/land (spikes) or stock/land stock/land might have spikes too, since dividing by land tricky at the changes, esp with lag

plotting_and_hector_code

acs

true

Setup relevant libraries and directory for saving figs and data

GCP data for NBP

Read in and understand GCP data

Check if GCP data and Dawn’s agree

Proceed with GCP data

compare hector NBP with GCP

LUC emissions received by Hector

Plot NBP data

Test Hector NBP constraints.

functions

Check effect of ocean param change

0-level test: Hector 245 1970-2100 constraint

0-level test: Hector 126 1970-2100 constraint

0-level test: Hector 585 1970-2100 constraint

More constraints

hector 245 1970-2021 constraint

GCP NBP constraint 1972-2021

245 runs magnified

hector 585 1970-2021 constraint

GCP 585 1972-2021

GCP 585 1972-2000

585 runs magnified

126 1970-2021

Constraints across scenarios

Conclusions

DISREGARD from here other tests

Ordered constraints 1970-2021 and T, CO2

Ordered constraints 1970-2050 and T, CO2

lower diffusivity 1970-2021

lower diffusivity 1970-2050

replicate NBP-T-C issue

sweep diffusivity vs gcp magnitude

test doing one year at a time (like DW work)